function matr = matrix_of_boundary_conditions(p, m, CONSTS)

    k0    = CONSTS.k0;
    a     = CONSTS.a;
    eta   = CONSTS.eta;
    eps_a = CONSTS.eps_a;

    %q = func_q_from_p(p, CONSTS);
    s  = func_s_from_p(p, CONSTS); 
    [q1 q2] = func_qk_from_p(p, CONSTS);
    Q1 = k0*a*q1;
    Q2 = k0*a*q2;
    S  = k0*a*s;
    
    bessel_m_q1 = besselmx(double('J'), m, Q1, 0); %besselj(m, Q1);
    bessel_m_q2 = besselmx(double('J'), m, Q2, 0); %besselj(m, Q2);
    bessel_mp1_q1 = besselmx(double('J'), m+1, Q1, 0); %besselj(m+1, Q1);
    bessel_mp1_q2 = besselmx(double('J'), m+1, Q2, 0); %besselj(m+1, Q2);

    bessel_K_m = besselmx(double('K'), m, S, 0); %besselk(m, S);
    bessel_K_mp1 = besselmx(double('K'), m+1, S, 0); %besselk(m+1, S);

    [n1 n2] = func_nk_from_p(p, CONSTS, false);
    %[J_m1 J_m2 J_m1_1 J_m2_1] = func_J_mk(q, m, CONSTS, false);
    [alfa_1 alfa_2 beta_1 beta_2] = func_alfa_beta_from_p(p, CONSTS, false);

    %K_m    = func_K_m(p, m, CONSTS);
    
    combo1 = bessel_mp1_q1 + alfa_1.*m.*(bessel_m_q1./Q1);
    combo2 = bessel_mp1_q1 - beta_1.*m.*(bessel_m_q1./Q1);
    combo3 = bessel_mp1_q2 + alfa_2.*m.*(bessel_m_q2./Q2);
    combo4 = bessel_mp1_q2 - beta_2.*m.*(bessel_m_q2./Q2);
    combo5 = bessel_K_mp1 - m.*(bessel_K_m./S);
    
    
%     matr = [Q1.*bessel_m_q1,              Q2.*bessel_m_q2,            0.0,                                     S.*bessel_K_m;  ...
%             n1.*Q1.*bessel_m_q1,          n2.*Q2.*bessel_m_q2,        (eta/eps_a)*S.*bessel_K_m,               0;     ...
%             J_m1.*Q1.*bessel_m_q1,        J_m2_1.*Q2.*bessel_m_q2,    (p./eps_a).*(m./(S.^2)).*S.*bessel_K_m,  -K_m.*S.*bessel_K_m;   ...
%             n1.*J_m1_1.*Q1.*bessel_m_q1,  n2.*J_m2.*Q2.*bessel_m_q2,  -K_m.*S.*bessel_K_m,                     (p.*(m./(S.^2))).*S.*bessel_K_m];

    matr = [Q1.*bessel_m_q1,      Q2.*bessel_m_q2,      0.0,                                     S.*bessel_K_m;  ...
            n1.*Q1.*bessel_m_q1,  n2.*Q2.*bessel_m_q2, (eta/eps_a)*S.*bessel_K_m,                0;     ...
            combo1,               combo3,              (p./eps_a).*(m./(S)).*bessel_K_m,        -combo5;   ...
            n1.*combo2,           n2.*combo4,          -combo5,                                 (p.*(m./(S))).*bessel_K_m];

end

% function y = func_K_m(p, m, CONSTS)
%     k0 = CONSTS.k0;
%     a = CONSTS.a;
%     s = func_s_from_p(p, CONSTS);
% 
%     y = (besselk(m+1, k0*a*s, 0)./((k0*a*s).*besselk(m, k0*a*s, 0))) - ...
%            (m./((k0*a*s).^2));
% end
